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The homogeneous dynamo effect is at the root of cosmic magnetic field generation. With only 
a very few exceptions, the numerical treatment of homogeneous dynamos is carried out in the 
framework of the differential equation approach. The present paper tries to facilitate the use of 
integral equations in dynamo research. Apart from the pedagogical value to illustrate dynamo 
action within the well-known picture of Biot-Savart's law, the integral equation approach has a 
number of practical advantages. The first advantage is its proven numerical robustness and stability. 
The second and perhaps most important advantage is its applicability to dynamos in arbitrary 
geometries. The third advantage is its intimate connection to inverse problems relevant not only for 
dynamos but also for technical applications of magnetohydrodynamics. The paper provides the first 
general formulation and application of the integral equation approach to time-dependent kinematic 
dynamos in finite domains. For the spherically symmetric o? dynamo model it is shown how the 
general formulation is reduced to a coupled system of two radial integral equations for the defining 
scalars of the poloidal and the toroidal field components. The integral equation formulation for 
spherical dynamos with general velocity fields is also derived. Two numerical examples, the a 2 
dynamo model with radially varying a, and the Bullard-Gellman model illustrate the equivalence 
of the approach with the usual differential equation method. 

PACS numbers: 



I. INTRODUCTION 



Cosmic magnetic fields, including the fields of planets, stars, and galaxies, result from the hydromagnetic dynamo 
effect P, Q . The last decades have seen tremendous progress in the analytical and numerical treatment of magnetic 
field generation in cosmic bodies. Recently, the hydromagnetic dynamo effect has been validated experimentally in 
large liquid sodium facilities in Riga and Karlsruhe 0, & IE • 

The usual way to treat hydromagnetic dynamos numerically is within the differential equation approach. Supposing 
the fluid velocity u to be given, the governing differential equation is the induction equation for the magnetic field B, 

^- = Vx(uxB) + AB, (1) 

Ot fJ,Q<T 

with hq and a denoting the permeability of the free space and the electrical conductivity of the fluid, respectively. 
Equation (JTJ follows directly from pre-Maxwell's equations and Ohm's law in moving conductors. Note that the 
magnetic field has to be divergence-free: 

V • B = . (2) 

In case of vanishing excitations of the magnetic field from outside the considered finite region, the boundary 
condition for the magnetic field reads 

B = 0(r~ 3 ) as r^oo. (3) 

The induction equation (1) is sufficient to treat kinematic dynamo models in which the back-reaction of the self- 
excited magnetic field on the flow is neglected. Such a simplification is justified during the initial phase of self-excitation 
when the magnetic field is weak. For stronger fields, one has to cope with dynamically consistent dynamo models which 
require the simultaneous solution of the induction equation for the magnetic field and the Navier-Stokes equation for 
the velocity. This saturation regime, however, will not be considered in the present paper. 
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The spherical geometry of many cosmic bodies such as planets and stars has simplified dynamo simulations in 
one important respect: for the spherical case the boundary conditions for the magnetic field can be re- formulated 
separately for every degree and order of the spherical harmonics, which makes any particular treatment of the magnetic 
fields in the exterior superfluous. 

When it comes to dynamos in other than spherical geometries, this pleasant situation changes. Then the correct 
handling of the non-local boundary conditions becomes non-trivial. Such a problem appears, e.g., in the numerical 
simulation of the recent dynamo experiments that are of cylindrical shape, but also in the simulations of galactic 
dynamos. There are some ways to circumvent this problem: one can use simplified local boundary conditions ("vertical 
field condition" ) H, E| , one can embed the actual dynamo body into a sphere with the region between the actual 
dynamo and the surface of the sphere virtually filled with a medium of lower electrical conductivity [Tol ITl| , or one 
can solve the Laplace equation in the exterior and fit the solution at the boundary to the solution in the interior |12| , 
which is, however, a tedious and time-consuming procedure. 

The integral equation approach that we will establish in the present paper is intended to change this unsatisfactory 
situation concerning the handling of boundary conditions. For the steady case, and for infinite domains of homogeneous 
conductivity, the integral equation approach was used in a few previous papers llfiL . The inclusion of 

boundaries, again for the steady case, was already delineated in the book of Roberts |l7| . In Roberts' own opinion 
C|l7|. p. 74), however, this formulation did "...not appear, in general, to be very useful." In [Tsj we have tried to 
put this pessimistic judgment into question. In particular, we have derived from the general theory a system of 
one-dimensional integral equations for a dynamo model with a spherically symmetric, isotropic helical turbulence 
parameter a in a finite sphere, and we have re-derived analytically the solution found by Krause and Steenbeck ppj 
for the special case of constant a. In |l9j | we have investigated the performance of numerical schemes to solve these 
integral equations for steady a 2 dynamos. 

It should be pointed out that a similar approach had been established earlier under the label "velocity-current- 
formulation" by Meir and Schmidt |2*lL 12^ . This approach had also the pronounced aim to circumvent the numerical 
treatment outside the region of interest. However, the numerical focus of this work laid more on steady, coupled MHD 
problems with small magnetic Reynolds number than on dynamo problems. 

In the present paper, we generalize the integral equation approach for steady dynamos in finite domains to the 
time-dependent case. In contrast to earlier statements on this matter [lil Il8l. fli^ that the use of the Green's function 
of the Hclmholtz equation will lead us to a non-linear eigenvalue equation, we present here a linear formulation of the 
eigenvalue problem. This makes the approach more attractive for numerical treatment. 

For the paradigmatic case of a time-dependent spherically symmetric, isotropic a 2 dynamo we derive, starting from 
the general formulation, a coupled system of two radial integral equations which is solved numerically. The equivalence 
with the results of an differential equation solver is made evident. 

We also derive the coupled system of radial integral equations for spherical dynamos with general velocity fields, 
and treat numerically the well-known Bullard-Gcllman model within the new approach. 



II. GENERAL FORMULATION 



In this section, the general form of the integral equation approach is developed for time-dependent dynamos acting 
in an electrically conducting, non-magnetic fluid which occupies a finite domain D with a boundary S, surrounded by 
non-conducting space. Under the assumption that the velocity or corresponding mean-field quantities are stationary, 
we generalize the basic idea and the methods from the steady case [3 to the time-dependent case. 

We start with pre-Maxwell's equations (the displacement current can be skipped in the quasistationary approxima- 
tion), 

V,E M = -«, (4) 

V-B(r,t) - 0, (5) 
VxB(r,t) = Mo j(r,i) , (6) 

where B(r, t) is the magnetic field, E(r, t) the electric field, r the position vector, t the time, and jiQ the magnetic 
permeability of the free space. The current density j(r, t) satisfies Ohm's law 

j(r,i)=tr(E(r,f)+F(r,t)) (7) 



inside the dynamo domain D, and it vanishes outside D. F(r, t) denotes the electromotive force (emf) u(r) x B(r,£), 
where u(r) is the velocity of the fluid motion, which we suppose to be stationary. In the framework of mean-field 
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electrodynamics (see, e.g., pj) B(r, t), j(r, t) and F(r, t) are split into mean fields and fluctuating fields. Then F(r, t) 
can be fixed to the form 

F(r, t) = u(r) x B(r, t) + o(r)B(r, t) - /3(r)V x B(r, t), (8) 

where u(r) and B(r, t) now denote the mean velocity and the mean magnetic field, respectively. The term a(r)B(r, t) 
describes the emf due to the non-mirrorsymmetric part of the turbulence (a-effect), and the term /?(r)V x B(r, t) 
describes another effect which can be interpreted as an conductivity decrease due to the turbulence (/3-effect). 

The equations given so far, together with initial conditions for B, define an initial value problem for B(r, f) which 
can be combined into the form of Eq. Together with the requirements that there are no surface currents on S 
and that B(r, t) vanishes at infinity, they allow to determine B(r, t) once the concrete dependence of F(r, t) on B(r, t) 
has been fixed. 

Under the condition that u, a and in Eq. JHJ only depend on the position r, the functions E(r, t), B(r, t), F(r, t) 
and j(r, t) may be written in the form: 

E(r, t) = E(r)e At , B(r, t) = B(r)e At , F(r, t) = F(r)e At , j(r, t) = j(r)e At , (9) 

where A is in general a complex constant. Then Maxwell's equations and Ohm's law transform into 

V x E = AB , (10) 
VB = 0, (11) 

V x B = /ioj , (12) 

j = <r(E + F) . (13) 

Hereafter, the symbols B, E, j and F denote the functions which only depend on the position r, but not on the time 
t. 

Since the magnetic field B is divergence-free, there exists a magnetic vector potential A(r) such that 

B = V x A, (14) 

which implies, together with Eq. (|10[) . that 

Vx(E + AA)=0. (15) 
Therefore, E + AA has to be irrotational, and we can write 

E + AA = -V<p. (16) 

Subsequently, Ohm's law gets the form 

j = c(F — AA - V(p). (17) 

In the following, we will derive a system of integral equations that is equivalent to the differential equation formu- 
lation. The starting point for the first integral equation is the application of Biot-Savart's law on Eq. 1)12(1. leading 
to 



= jo , ^l^ZU dV > . (18 ) 
4tt J d |r-r'| 3 

It is well known that the curl operator in Eq. I|12[> is a left inverse of the Biot-Savart operator when the current j 
is divergence-free and tangent to the boundary of the domain |2^|. Both conditions are indeed fulfilled in our case. 
Inserting Eq. I|17fl into Eq. Q18|). using Eq. I|14|) and Gauss's theorem, we get the first integral equation 

= W r F(rQ x (r - rQ ^ _ ^A t B(rQ ^ 
4tt J d |r-r'| 3 4tt J d |r-r'| 

Mo^A f , A(s') fi o- f / /\ / /\ r 



-in 



J ■«) x m*-!£j ,(On(0 x ^ iS ', (19) 



with n(s') denoting the outward directed unit vector at the boundary point s' and dS' denoting an area element at 
this point. For some purposes it might be useful to express the first volume integral in Eq. i|19|) i n the form 

^ dv = / - 1 „(o *m, ne . (20) 

|r-r'| 3 J D |r — r'| J s |r-s'| 
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For the steady case (A = 0) Eq. 1(19(1 reduces to the form given in [T^, with one volume integral over B and one 
boundary integral over tp. The latter one would vanish only in the case that the dynamo region is extended to infinity. 
The time-dependence introduces now a new volume integral over B, but also one boundary integral over the vector 
potential A which cannot be reduced to a simple expression in B. Before we focus on this point, let us first derive 
the integral equation for tp. 

From Eq. 1(17(1 and the demand that the current has to be divergence-free, V • j = 0, we get a Poisson equation for 

tp: 

Aip = V • (F — AA) . (21) 

Assuming vacuum boundary condition which means that the current must not leave the domain D, we obtain from 
Green's theorem the following boundary integral equation for tp: 

If F(r') • (r - r') A f V r , ■ A(r') 

P ip{T) = — / j jr. dV + — — — dV 

47r J |r-r'| 3 47r J |r — r'| 

D D 

s s 11 

where p = 1 for points r inside D, p = 1/2 for points r = s on S and p = for points r outside D. Again, another 
expression of the first volume integral might be useful : 

|r-r'| 3 J D |r-r'| J s |r - s'| 

If we use the Coulomb gauge for the vector potential, V • A = 0, we see that the second volume integral in Eq. I|22|) 
vanishes. Of course, one is free to use other gauges than the Coulomb gauge. One advantage of this gauge, in view 
of later numerical applications, is that it leads to a formulation in which the vector potential is only needed at the 
boundary of the domain. 

For the steady case, Eqs. (|19|) and l|22l) with A = are sufficient to determine the magnetic field B. But for 
the time-dependent case presently under consideration, we have to introduce the vector potential A, at least at the 
boundary, for completely formulating the problem. Necessarily we have to establish another relation for A in order to 
make the problem solvable. From Eq. (|14fl and Helmholtz's theorem ([2^, p. 53) we can express the vector potential 
at the boundary in one of the two forms 

aw = -L / l^War 

iwJ D |r - r'| 

4ttJ d |s-r'| 3 AttJs \s - s'\ 

The integral equations 1)19(1 . 1)22(1 and 1(24(1 provide another complete formulation of the problem for B. The main 
advantage of this formulation is that one can avoid any treatment of fields in the exteriour of D. The boundary 
conditions are being fulfilled by solving the boundary integral equations for ip and A. 

In the next section we will sketch how the equations for A and tp at the boundary can be absorbed into a single 
integral equation for B. 



III. GENERAL NUMERICAL ASPECTS 



In this section we delineate the general framework for the numerical solution of the coupled equations l(19|) , 122|) 
and l|24(l. Let us assume certain spatial discretizations of the magnetic field in the volume of the dynamo, and of the 
electric potential and the vector potential at the boundary. Then, Eq. I(19|) may be rewritten in the form 

B t = L lk B k + \M lk B k + XP in A n + N a tpi , (25) 

where the Bi denote the degrees of freedom of the magnetic field in the volume of the dynamo, while tpi and A n denote 
the degrees of freedom of the electric potential and the vector potential at the boundary. Here and in the following 
we use Einstein's summation convention, and we reserve the indices i, k for magnetic field degrees of freedom in the 
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FIG. 1: The general scheme for the numerical realization of the integral equation approach. Left: The connection between the 
magnetic fields in the volume and the vector potential and the electric potential on the boundary. Right: Reduction to the 
pure connection between magnetic fields in the volume. The boundary on the l.h.s. is indicated to have a finer discretization 
then the volume. On the r.h.s. the influence of the boundary has been incorporated implicitly. For details, cp the formulae in 
the text. 



volume of the fluid, whereas the indices I, m, n are reserved for the vector potential and electric potential degrees of 
freedom at the boundary of the fluid. 

For any given dynamo source, any shape of the dynamo domain, and any concrete form of the discretization, the 
matrices L, M, N and P in Eq. I|25|) can easily be derived from Eq. Q19[l. It is worthwhile to note that only L 
depends on the dynamo source (u or a), whereas M, N and P depend only on the geometry of the dynamo domain 
and the discretization details. 

Similarly, the discretization of the boundary integral equation 1221) (for the case that r is on the surface S) leads to 

0.5 ipi + E hn tp m = Hi k B k + XDi n A n , (26) 

or 

Gimfm = HikBk + XDi n A n , (27) 

where G/ m = 0.5 5i m + Ei m . Again one should note that only H depends on the dynamo source, whereas G and D 
depend only on the geometry of the dynamo domain and the discretization details. The discretization of Eq. I|24|) 
gives 

A„ = QnkBk , (28) 

with Q depending solely on the geometry. The influence of the magnetic field degrees of freedom on the vector 
potential and the electric potential degrees of freedom on the boundary, and vice versa, is illustrated on the l.h.s. of 
Fig. CD 

Substituting Eq. into Eq. J23 yields 

fm = {G 1 )mlHikBk + \{G~ 1 ),niDi n Q n kBk- (29) 
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However, for the inversion of the matrix G some care is needed. Basically, G is a singular matrix, reflecting the 
fact that the electric potential is determined only up to a constant. Accidcntly, it may happen that this singularity 
is weakened by inaccuracies due to the discretization. Nevertheless, one should be careful with the inversion. A 
convenient method to deal with the inversion is the application of the deflation method p4l l25j|. In the following, we 
will simply assume that the inverse of G has been found in an appropriate manner. 
After inserting Eqs. J2HJ and Eq. (|2"5|) is transformed to 

Bi = LikBk + XMikBk + XPinQnkBk + N im (G~ 1 ) m iHikBk + \Nu(G~ 1 )i m D mn Q n kBk- (30) 

Evidently, the electric potential and the magnetic vector potential at the boundary served only as auxiliary quantities 
in order to ensure the right boundary conditions. From the numerical viewpoint it is important to take notice of the 
following: for the accurate solution of the boundary integral equation l|29|) for a dynamo in a given domain it may be 
advisable to use a fine discretization of the boundary with a large number of grid points. Hence, the corresponding 
inversion of the matrix G might be numerically expensive. However, for a given geometry this inversion is needed 
only once. Finally, after carrying out the matrix multiplications N ■ G 1 • H and N • G 1 ■ D • Q in Eq. (|30(l one 
ends up with a matrix of the order (NB,NB) where NB denotes the total number of all magnetic field degrees of 
freedom (see the r.h.s. of Fig. Q). 

Now, Eq. H3Uf) can be rewritten in the following form: 

(6 ik - L ik - Ni m {G- l ) m iHi k )B k = X(M ik + P in Qnk + N l i(G~ 1 )i m D mn Q n k)Bk- (31) 

This is a generalized linear matrix eigenvalue problem in which only the magnetic field components remain. The 
numerical solution of the arising linear generalized eigenvalue equation 1)3 l|l yields the eigenvalues A, comprising as 
the real part the growth rate and as the imaginary part the frequency of the dynamo. 



IV. APPLICATION TO TIME-DEPENDENT a 2 DYNAMOS - BASICS 



In this section, we exemplify the general integral equation approach by applying it to a simple mean-field dynamo 
model with a spherically symmetric, isotropic helical turbulence parameter a. In contrast to the original model with 
constant a [H l20j|. whose advantage is the possibility of an analytical treatment, we allow here a to vary with the 
radial coordinate r. Starting from equations (|19|1 . (|22H and H24|) we will derive two coupled radial integral equations 
which will be used in the next section for numerical treatment. Note that for the steady case the corresponding 
derivation had been given in [Isj . 

Basically, this section is intended as an exercise to show how the basic integral equation approach can be reduced 
in a particular case with high symmetry of the dynamo source. A slightly simpler derivation of the two radial integral 
equations, which starts from the radial differential equations and uses a Green's function method, is given in Appendix 
A. If the reader is only interested in the final form of the radial integral equations he can skip the following derivations 
and refer directly to Eqs. (|T7jl and (|(77|) . 



A. Preliminaries 



As usual in dynamo theory, we split the divergence-free magnetic field B into a poloidal and a toroidal part, denoted 
by Bp and Bj. Since we use the Coulomb gauge, V • A = 0, an equivalent decomposition can also be applied to the 
vector potential, A = Ap + Ap. We represent these fields by the defining scalars S, T, S A , T A according to 

B P = VxVxf~rV B T = Vx^rj , (32) 
A P = V x V x (— rV A T = V x H— r\ . (33) 
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We introduce spherical coordinates r, 9, <p and denote the radius vector by r. The defining scalars and the electric 
potential are expanded in series of spherical harmonics Yi m (9, (f>), 



S(r,6,<f> 
T(r,6,<j> 
S A (r,0,cl > 
T A (r,6,(f> 
<p{r,6,<p 



mi"] Yh 



^ sim(r)Yi 

I, m 
Lin 

y^(/?; m (r)Yj m (6>, 



(34) 



Later S , T , ip will be needed only at the boundary r — R. For the spherical harmonics Y/ 



lm \ 



the definition 



121 + 1 (/-m)! 
47r (l + m)\ 



P lm {cos9)e lm * 



(35) 



is employed, with P\ m denoting associated Legendre Polynomials. The summations in l|34(l are over all degrees I and 
orders m satisfying I > and |m| < terms with 1 = 0, however, are without interest in the following. Since S, T, 
S A , T A , and are real we have s;_ m = s* lm and analogous relations for i; m , s^, t A m and v?im- The definition 135|) 
implies the following orthogonality relation for the Y\ m (9, <fi): 



2tt 



A useful relation is 



where the operator is defined by 



sin0d0y,;T m ,(M)*J m (M) = 5 w S n 



ft Yi m — —1(1 + l)Yi m , 



1 d 



Of 



09 



1 d 2 f 



sm 



From Eqs. i|32|) , l|33|l . and l|34|) we obtain, with the help of (|37|) . the components of B 
Br(r,0,4>) 



^ - ^ SZm(r)Yi m (fl, 0) 



E 
E 

Lm 



ti m (r) dYi m (6,(f>) 1 dsi m (r) dYi m (9,< 



r sin r dr 90 

Mr)dY im (0,0) , 1 d.s im (r) 9^ m (0, 



9(9 



r sin dr 



and of A 



MrAti = J2 l -^2^ s t(r)Yim(9,t) 



A 9 (r,9,q>) = Y, 



l,m 



, ldsf m (r)dY lm (0,< 



A^(r,9,q>) = E - 



Lm 



r sin # 90 r dr <9# 

tL<r)3Yi m (9,<t>) , 1 ds^(r)SH, 



r sin dr 



(36) 



(37) 



(38) 



(39) 



(40) 
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Finally, we recall the expression for the inverse distance between two points r and r', 

u— ^ = 4 "E E —^YLie'^YUe^) , (4i) 



1=0 m=-l 



where r > denotes the larger of the values r and r', and r< the smaller one. 

Equipped with these preliminaries, we will derive two coupled integral equations for the functions si m (r) and 
U m (r). To meet this goal, we will also need expressions for s^^r) and tp(r) which are necessary, however, only at 
some intermediate steps. 

B. Radial integral equation for s; m (r) 

Taking the scalar product of both sides of Eq. I|19|l with the unity vector e r we obtain 

R , , W f Mr')B(r') - AA(r')) x (r - r') 
B ( r ) ' e r = — ^ ^—^- 3 e r dV 

' I s ^ (s>(s,) x ]7— ^ ■ e " dS ' 

■ e r > — dV . (42) 



47T 

Ho<r f V w x (a(r')B(r') - AA(r')) 



s 

r' 

'■ — • fi_y 

47T |r — r 

In the derivation of the last step in Eq. (|42|l we have expressed e,. under the integrals by (r — r')/r + (r'/r)e r /, and 
we have used the fact that the triple product vanishes when n(s') and e r < coincide for r' = s'. 
By virtue of Eq. (jTi|) , Eq. becomes 

= ^ r V, x HrW)) ✓ /" B(r^ 

W 4vr 7 D |r-r'| r 4tt 7 d |r-r'| r V ' 

Noting that in Eq. |g2J) 

V r / x (a(r')B(r')) = -B(r') x V r /a(r') + a(r')V r , x B(r') , (44) 

we see that the scalar product of the first term on the right hand side with e r / vanishes since the gradient of a(r') 
points in r'-direction, too. From Eqs. 13 7|) and 139fl we obtain 

(V r / x B(r')) ■ e r , = ]T l ^±^-t Vm , [r')Y Vm , (<?', 0') . (45) 
Taking Eqs. g2J>> (01} and 03]) together we find 



\ - + 1) ,s v ,x A*oo- / , \ - l'{V + 1) 4 , 
2^ -^2— ^ m (r)F im (^,^) = — y a(r )^ ^ 2 frw(r , r r , 



M0 aA y ^ V^^^y^ft*^ . (46) 



47rr 



After expressing the inverse distance according to Eq. (|4*T1) , integrating on the right-hand side of Eq. (|4^|l over the 
primed angles, multiplying then both sides of Eq. (|46|l with Y* m (9, 4>) and integrating over the non-primed angles we 
obtain the first integral equation of our problem in the form 

s im(r) = ^T[l ri —ra{r')t lm {r')dr'+ /" ^ a(r') t lm (r') dr' 



21 + 1 



r' 1 



r r ll+l t _ rR r l+l 



si m {r')dr' - A / — Slm ( r ')dr'} . (47) 
r 
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C. The electric potential at the boundary 

For the determination of the electric potential at the boundary it is convenient to start from Eq. (|22H for points r 
outside D. As for the last boundary integral in Eq. IL'l'l) we have 



±- / ^(s')n(s') • r S ', dS' = ±- I V tf)-!L—L- dS' 
4ttJ ^ v |r-s'| 3 47ri ^ ; 9s'|r-s'| 



s s 



J J2wrn(R)Y lm (6',<p')J2 y^^^ Y i* m >(0'A')Yi> m '{e,ct>) dS' (48) 



and thus 



lim j- [ ^(s')n(s') ■ dS> = J2 ^-rV>lm(R)Y lm ( 

q 1 1 lm 



(49) 



For the evaluation of the volume integrals in Eq. I|22(l we can make the second one vanish by means of the Coulomb 
gauge V • A. For the first one, we use the alternative formulation Eq. I(23j) . Taking B r from Eq. Ij39j) ■ we find 

V r , • Mr')B(r')) /■ da(r') ^ i'(Z' + 1) 1 

1 1 Z'm' 1 1 

and thus 



1 / V r ,-( a (r')B(r')) ^ v^ jq + l) v m f ^ da(r>) 

HtoJ jF=?j ^ = L^T^'H i^+T ~dP~ 

j-) lm 



si m {r') dr . (51) 



Analogously, we obtain for the remaining boundary integrals in Eq. 11' 2 1) 

,. If,,, a(s')B(s') - AA(s') ^1(1 + 1)1. ._. x A 

(52) 

Evaluating now Eq. for r — > s with the help of Eqs. |@3J|, (f3T~f> and we find 

f R r' 1 da(r') I + 1 

<p ha (R) = -(l + l)J ^-^s lm {r')dr' + — (a{R)s lm {R)-\s? m {R)). (53) 

This expression for <p(R) will later be needed in the integral equation for ti m (r). 

D. The vector potential at the boundary 
From Eq. H24[) and the fact that the curl of the magnetic field vanishes outside of D, we obtain 

v -; xB( , r V = / V B fV, (54) 



'D I 1 _ 1 I JD+D' I 1 - 1 I 

where D' denotes the outside region of D. Note that 

V r » x = V,.^- x B(r') + ^^V r , x B(r'). (55) 

|r — r'| |r — r'| |r — r '| 

The application of this equation on the right side of the Eq. 154(1 leads to 

^C^ dV' =[ V r , x - [ V,* x B(r')^'. (56) 



D+D' I 1 — 1 I JD+D 
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Applying Gauss' theorem, we have 

I V^xB(O r = I x W ^ I r-r_ x 

Jd F-r'| 7 Soo |r — r'| 7 D+D , |r - r'| 3 

Equation Q allows to conclude that the surficial integration on the r.h.s. of this equation vanishes. Therefore, after 
taking the scalar product of both sides of the equation l|57|l with the unit vector e r , we obtain 

I Vr ,' X B( , r,) • e r dV = f V r ,^- x B(r') • e r ,-dV. (58) 
Jd |r-r'| J D+D , |r-r'| r 

In the derivation of this equation, the relation e r = (r — r')/r + (r 1 /r)e r > has been used again. Applying Eq. I|55|) 
again we have 

V r / x B(r') , . f x B(r') r' , T . /" „ B(r') r' „ rl 

IZ-ZFi ■ *rdV = / ■ e r / -dV' - / V r / x 7—^-4 • e r / -df' (59) 



The second term on the r.h.s. of Eq. (|59J) can be shown to vanish, so we conclude that 



A(r).e r = -L m^l.e/- d V. (60) 
47r In r — r r 



Combining Eqs. (|4T?|) . and (|4l>)) . we obtain 



*Lir) = 5r - T ( / -j-^CrOdr' + / - r i Jm (r')iir'). (61) 



2/ + 1 Vo 

Particularly, when r — R in the above equation, we have 



j pR r ,l+i 



*L(R) = J o -^t lm (r')dr'. (62) 



As we will see, there is no need to calculate the corresponding expression for tf m (r) . 



E. Radial integral equation for tt m (r) 

We take now the curl of both sides of Eq. (|18fl thus obtaining 

V r x B r = — V r x V r x / , — ; —dV 

4tt J d |r — r'| 

-V P x / ^(s')n(s') x dS> ]. (63) 

Js |r-s'| J 

Considering first the case r < R we take on both sides of Eq. (|63[) the scalar product with e r . We note that 
er • (Vr x V, x [ W)-AAy) r 



e,((V r V,-A r )/ ;W)-^ r) 



d f V r / ■ (a(r')B(r')) rfy , _ 3_ a(s')B(s') - AA(s') ^ 



<9r /d |r — r'| dr 

s 

+4ir(a(r)B r (r) - \A r (r)) , (64) 

where we have used the identity A r |r — r'| _1 = — 47r<5(r — r') and the Coulomb gauge V • A = 0. The two integrals 
on the third line of Eq. I|64|l were already treated in subsection C. Concerning the boundary integral in Eq. I|63|) over 
the electric potential we note that 

Br ■ ( V ' X ( n(S ' } X Tr^p)) = -£^vh\ ■ (65) 
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Putting everything together we obtain 



(V r x B(r)) • e r - fx a(a(r)B r (r) - XA r (r)) + 

47T 



drJ D dr' rl j |r-r'| V 



8 f - XA r (?) ^ + f # 1 ^ 



s 



s 



Representing now the left-hand side according to Eq. 1(35}, applying Eq. I|41l) and integrating both sides over the 
angles we obtain 

, / \ r i \ i \ ^ a , \ I + 1 f r da(r') . ,.r' , . 

ti m {r) = Hov[a(r)s lm (r) - Xs lm (r) - ^p-j — s im (r )- j -dr 

I f R da(r') . .r l+1 J , I r l+1 . /m , m 
-A ^ m (i?)) - ^^(i?)^] . (66) 



Substituting Eqs. J53J| and into Eq. JBBJ leads to 



l + l r da(r') . ,,r", , 



Um(r) = fM)a[a(r)si m (r) ~ 2 i + \ J dr' Slm ^ r ^~p: dr ' 

I f R da(r>) r l +^ l + l r l + l f R ,ida(/) 



A r'+! /■« , „ . , A rr ,,+1 



ti m {r')dr' 



\ rR l+l l+l 
~ 21+1 J — tl ™^ dr ' - rt+T a ( R > l ™(R)} ■ (67) 

This is the second radial integral equation for our problem. Therefore, the spherically symmetric a 2 dynamo model 
is reduced to the two coupled integral equations l|47(l and l|67|l . 

F. Connection with the differential equation approach 

Notwithstanding the fact that the differential and the integral equation approach are equivalent in a general sense 
it might be instructive to show this equivalence for our special problem. Differentiating equations (|47|) and l|67|l two 
times with respect to the radial component the following relations can be obtained: 

1 d 2 si m 1(1 + 1) 

Xsim = — — = ^ — sim\ + a(r)ti m , (68) 

[1q(7 dr^ r z 

. , 1 f d 2 t lm 1(1 + 1) d ds im 1(1 + 1) 

XU m = — — - — ti m \ - - a r - — ) H 5 — a(r)s; m , (69 

finer dr 1 r z dr dr r z 

t lm (R) = R^^l\ r=R + l Slm (R) = 0. (70) 
dr 

As expected, these are the differential equations for the considered problem of radially varying a for the time-dependent 
case |2(| [27) • Note again that from the radial differential equation we can also derive the radial integral equations by 
means of a Green's function technique (see Appendix A). 
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V. APPLICATION TO TIME-DEPENDENT a 2 DYNAMOS - NUMERICS 

A. Discretization 

In the previous section, and in Appendix A, we have derived the radial integral equations (j47(l and (|67|l governing the 
time-dependent dynamo problem. In this subsection, we develop a numerical method to solve this integral equations 
system. 

Let us introduce the following definitions: 

x = r/R,Ca(x) = R 2 (ioaa(x) , A/ = R ^aXi, (71) 
where C is the magnitude of the function R 2 fiQaa. In addition, we introduce the notations 

G °^=\-^SOx>2' (72) 

Gt(a; '" o) -\^(^ +1 -^)4 +1 ^>-o, (73) 

ri T T „\ - J 2/+l X A "t- 2 /+l x x 



2Z+1 ~ 2Z+1 ^0' — U' 



Then, Eqs. (g7|) and ((HZI) obtain the form 



fi _ pi 

si m (x) = -C I G s (x,x )a(x )ti m (xo)dxo + Xi G s (x 1 x )si m (x )dx , (75) 



Jo 

ti m {x) = Ca(x)si m (x) -Cx l+1 a(l.Q)s lm (l.Q) + C [ da j x ° } Slm ( Xq ) Q t (x ,x ) rfx 

Jo "^o 

+A; / G t (a;, x )ti m (x )dx - C f x~ l x L da j- x °^ Slm (x a )dx . (76) 
Jo Jo " x o 

Choosing N equidistant grid points xi — iAx with Ax — l/N and approximating the integrals by the extended 
trapezoidal rule according to 



,1 n 

/ f(x)dx* V-(/(^_ 1 ) + /(ar i ))A^ 



(77) 



we obtain 



2V W-l 

A; si m (xj)G s (xj, Xj)Axcj = si m (xi) + C a^^iim^^G^i,-, a;j)cjAx, (78) 



AT-l 



A/ Gt(xi,Xj)ti m (Xj)AxCj = ti m (xi) - Ca(i,)s lm (3;,) + Ciz;[ +1 a:(1.0)s; m (1.0) 



„-^da(xj) 
2^ — ^ — s; m (x :) )Gt(xj,Xj)Aa;c : ,- 

+C^xr'4^^1 Sim (x J )Axc J , (79) 
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TABLE I: Comparison of the calculated growth rates and the analytic ones for the free decay case a(x) = 0. The degree of 
the spherical harmonics is I = 1. n = 1...4 correspond to to modes with increasing radial wavenumber n. The last row shows 
the analytic results. The other rows express the numerical results obtained by the integral equation approach for different grid 
numbers N . 



N 


n=l 


n=2 


n=3 


n=4 


8 


-9.79494 


-37.71179 


-79.61435 


-129.36969 


16 


-9.85079 


-39.02608 


-86.41265 


-150.20904 


32 


-9.86485 


-39.36475 


-88.21583 


-155.94894 


64 


-9.86844 


-39.44979 


-88.67356 


-157.42004 


128 


-9.86933 


-39.47085 


-88.78596 


-157.78870 


analytic 


-9.86960 


-39.47842 


-88.82644 


-157.91367 



where cm = 0.5, and Ci — 1.0 for i = 1, 2, • • • , N — 1. Equations l|78|l and l|79|l can be written in the matrix form: 

A;VX = WX, (80) 

with 

X = (si m (xi),Si m (x 2 ), - ■ ■ ,Sl m (x N ),t lm (xi),ti m (x 2 ),--- ,tlm(xN~l)) T , 

Vi.j = G s (xi,Xj)Axcj, = 1,2, • • • ,N) 

V l+N .j = Q,(*=l,2,.-.,AT-l,j = l,2,.-.,Af), 

Vi, j+N = Q,(*'=l,2,---,iV,j = l,2,---,iV-l) 

Vi + N, 3 +N = G t (xi,Xj)cjAx, (i = 1,2, ••• ,JV - 1, j = 1,2, ••• ,JV — 1), 

Wij = 5{i,j),(i,j = 1,2,- ■■ ,N) 

W i}j+N = G s {x l ,x j )t\xc j a{x j ), (i = 1,2, ■ • ■ ,N, j = 1,2, • • ■ ,N - 1) 

Wi+Nj+N = S(i,j),(i,j = 1,2,- •• ,N-l), 

Wi+Nj - -C^^-Gtixux^Axcj ~Ca(x l )8(i,j) + Cx l + 1 a{l.Q)8{j.N) + 

Cxr l x]^^-Ax Cj H(i- j),(i= 1,2, - ■ ■ , JV - 1, j = 1,2,- • • ,N), (81) 
where the functions S(i,j) and H(i — j) are defined as 

S(i,j) = 

and 

H{i-j) 



1,« =j, 



1,*>J, 
0, i < j. 

Note that Eq. I|80fl is a linear generalized eigenvalue problem. By multiplying both sides of Eq. I|80|l by the inverse 
of the matrix V, we can convert it to the following standard eigenvalue problem: 

V _1 WX = A ; X. (82) 

This eigenvalue problem can be solved by standard numerical routines. First, the matrix V _1 W is reduced to the 
Hessenberg form, then the QR algorithm can be employed to obtain the eigenvalue A;. 



B. Numerical results 



In this subsection, we illustrate the numerical performance of the integral equation approach formulated in this 
paper by a few examples for the functions a(x). 

Let us start with the case a(x) = 0, which corresponds to a pure field decay within a conducting sphere. For this 
case, the eigenvalues A; are known from quasi-analytic calculations Q. In Table {1} (for / = 1) and Table {nj (for 
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TABLE II: Same as Table □ but for I = 2. 



N 


n=l 


n=2 


n=3 


n=4 


8 


-19.87668 


-55.87647 


-103.28617 


-155.45810 


16 


-20.11100 


-58.68925 


-114.71107 


-186.07856 


32 


-20.17073 


-59.42935 


-117.83333 


-194.82947 


64 


-20.18561 


-59.61664 


-118.63189 


-197.09552 


128 


-20.18911 


-59.66311 


-118.83366 


-197.66728 


analytic 


-20.19073 


-59.67951 


-118.89986 


-197.85780 
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FIG. 2: Relative error of the numerically determined eigenvalues for the free decay case, C = 0, I = 1, n = 1...4. The 
convergence behaves like iV -2 . 



I = 2) we have listed the numerical results of the integral equation solver for the eigenmodes with increasing radial 
wavenumbers n — 1...4, together with the analytical results. From these tables it becomes obvious that even for a 
few grid numbers robust results can be obtained. In Figs. [3 and we have plotted the relative errors of the results in 
dependence on the used number of grid points N. As it is typical for integral equations of that kind 0], the relative 
error decreases like N~ 2 . 

Another quasi-analytic result exists for the steady case Q: for I = 1 and C = 4.4934095 or I = 2 and C — 5.7634593 
we know that the first eigenvalues have to be zero. Table shows the results of the integral equation solver, again 
for I = 1 and / = 2, but only for n = 1. The convergence of the results, which is again ~ N~ 2 , is depicted in Fig. 0] 



15 
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FIG. 3: Same as Fig. El but for I = 2. 
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FIG. 4: Convergence of the first eigenvalue (n = 1) to zero in the cases: C = 4.4934095, / = 1 and C = 5.7634593, / = 2. The 
normalization of the error is done with respect to the growth rates at C = which are Ai — —9.8696 and Ai — —20.1907. 



Now, we turn to a more complicated case. It corresponds to the profile a(x) = C(— 21.46 + 426.41a; 2 — 806.73a; 3 + 
392.28a; 4 ). The choice of this somewhat strange function is motivated by the fact that it is an example of a proper 
oscillatory a 2 dynamo j^]. In Fig. [S]we show the results of the integral equation solver for the case I = 1 and n — 1...7. 
We see that the spectral dependence on C (Fig. |SJ is very complex, with merging and splitting points of neighboring 
branches at which non-oscillatory solutions turn into oscillatory solutions and vice-versa. The computation was done 
with a grid number of N = 128, and the result is basically identical with that of a sophisticated differential equation 
solver ,27j. Hence, Fig. [S] might serve as a striking example that the integral equation approach works satisfactorily 
also in case that complex eigenvalues appear. 
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TABLE III: Convergence of the first eigenvalue (n = 1) to zero in the cases: = 4.4934095, I = 1 and = 5.7634593, 2 = 2. The 
normalization of the errors is done with the growth rates at C = 0, which are Ai(C = 0) = —9.8696 and \2{C = 0) = —20.1907 



1 N=8 N=16 N=32 N=64 N=128 

1 0.19836 0.04971 0.01249 0.00355 0.00052 

2 0.58527 0.14903 0.03740 0.01010 0.00072 




n=2 /' \ 


n=7 

n=4 n =3 




n=1 V V 


I f~~ ' 

n=6* 





0.5 1 1.5 2 



FIG. 5: Special case a(x) = C(-21.46 + 426.41x 2 - 806.73a: 3 + 392.28a; 4 . Growth rates and frequencies for the eigenfunctions 
with (2 = 1, n — 1...7). The number of grid points was iV = 128. Note the merging and splitting points of the spectrum, 
indicating transitions from non-oscillatory to oscillatory modes and vice versa. 
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VI. GENERAL VELOCITY FIELDS IN SPHERICAL GEOMETRY 

A. Basics 

In this section the Green's function method from Appendix A is applied to convert the induction equation for 
general velocity fields in a unit sphere to the integral equations system. In dimensionless form, Eq. (|TJ can be 
rewritten as follows 

/)B 

— = R m V x (u x B) + V 2 B, (83) 
at 

where R m is the magnetic Reynolds number, u and B may be expanded into the following series (Bullard and 
Gellman, 

u = ^(t a +s a ), (84) 

a 

B = ^(T^ + S^), (85) 



where 

t a = V x [e r t a (r,t)Y a (e,<p)], (86) 
s a = Vx Vx [e r s a (r,t)Y a (d,ip)], (87) 

etc. From here on, Y a denote the (2a + 1) surface harmonics P ama (9)sin(m a ip), P a m a (0)cos(m a ip)(m a — 0, • • • , a), 
where P am , a is a Legendre function with Neumann normalization, P a o the Legendre polynomial. Similarly s a 
is an abbreviation for s™°, etc. The summations in ()84|l and l|85|) are over cosine and sine contributions, 
a = 1, 2, 3, ■ ■ ■ : m n = 0, • • • , a; and similarly for (3,mp. Substituting Eqs. (|84|l and l|85|l into Eq. Q83JI . Bullard 
and Gellman (|3lJ) derived the spectral form of (|83|l as 

^ dSl 1 -^^S^^Y,^ S ^) + MS,) + ( Sa S p S 1 )i (88) 



a, (3 



dr 2 dt r 

~ ^ - = ~T Et^W) + (taS T 7 ) + (s a TpT 7 ) + [sctSpTj], (89) 



where 





— c\t a Sp, 










= c 2 s a Tp, 










dSp , 

= C3S a — h 

Or 


ds a 


j 






= c 5 t a Tp, 








(t a SpTy) 


dSa 

= c 6 t a — h 

Or 


dt a 

c?(^ 

or 


2/ 

— )Sp, 
r 






dTp 

= c 8 s a — h 

or 


ds a 


— )Tp + c 9 
r 


ds arr 
-dr- Tf) - 


(SaSpT^) 


d 2 S 


ds a 


dSp s a 


9Sp ^ 




+ cn a 
Or 


-5 h C12 — 

or r 


dr 



(90) 



<9 2 s Q 2 9s a 
c i3(^--^)^- 



The constants Cj in Eq. l|9"U|l are defined as follows: 
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Cl 



(-2 



W + i) 



iV 



La(a + 1) 



Ar- 



cs - -^a(a + l)(a(a + 1)- 0(0 + 1)- 7(7 +1)), 
c 4 = -^-0(0 + l)(a(a + 1)- 0(0 + 1)+ 7(7 + 1)), 



C5 



L 7 ( 7 +l) 



iV 7 



c 6 = - JL (0(0 + l)( a (a+l)- 0(0+1) 

+7(7 + 1)) + 7(7 + + 1) + 0(0 + 1) " 7(7 + 1))), 

c 7 = --|-0(0+l)(a( a +l)- 0(0+1) +7(7+1)), (91) 

c 8 - —a(a + l)(-a(a + l)+ 0(0 + 1) +7(7 + 1)), 

eg = — 7(7 + 1) K" + l)+ 0(0 + 1) -7(7 + 1)), 

cio = ^7-"( a + 1), 

en = ^-(a(a + l) + 0(0 + 1) -7(7 + 1)), 
2L , 

C12 = a+a + 1), 

C13 = ^0(0+D. 

In Eq. (J5TJ we have used the expressions for the Adams-Gaunt and Elsasser integrals 



K = I / YaYpYysinedOdcp, 
Jo Jo 



Jo 

27T />7T 



JO 



,8Y 8Y 1 dYpdYy 



L= y a (-^-^-^-^)d6d^ (92) 



89 dip dip 89 



and the normalization factor 



AT I 27+1 (7-m)!' 7" , 

7 - 1 4^ 7 ( 7 +l) Q 

L 27+1 ' 111 u - 

By the same Green's functions method as in Appendix A, and using the definitions for G s (r, ro), Gt(r, tq) there, 
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the differential equations system l|88|l and (|89|l can be converted to the following integral equations system: 

S -y( r ) = y^, R m I Gs ^ r ' r °" > [cit a (ro)Sd(r ) + c 2 s a (ro)Tp(r a ) + c 3 s a (r Q ) ^ 
1$ Jo r " dr ° 

ds a {ro) 



+c, 3 (%^--^)S 3 (ro)]dr„ + A / <7, ( ,-. n, )7i ( ,„ ulr n . ,94) 



+C4 — - Sf3{r )}dr + X G s (r,ro)S 1 (ro)dro, (93) 

u- r o Jo 

rp { \ d f 1 G t {r,r ) ,dS {r o ) f dt a (r ) 2t a (r ) 
T i( r ) = ) R™ § Ma( r o) r /3( r o) + c 6 £a( r o) — j 1- c T ( — )5/3(r ), 

+c 8 s a (r ) — h c 8 ( — j s Q (T-o))T^(r ) + c 9 — I>(r ) 

dr dr r dr 

d 2 Sp{r ) , ds a (r a ) dSp(r ) s a (r ) dSa(r ) 

+ CloS Q (ro) j-n 1" ClX 3 3 h C12 3 

drg dr dr r dr 

^ 2 s a (r ) 2 ris a (r c ^ ^ 
dr^ r dr Q 

Strictly speaking, this is an integro-differential equation system which could be used for numerical analysis. If one 
would insist on having a pure integral equation system one could employ integration by parts in order to obtain: 

S 7 {r) = Y^Rmlj G s {r,r )F 1 (ro)Sp(r )dr + f G s (r,r )F 2 (r )T^r )dr 

a/3 J ° J ° 

+ I dG i r,r0) F 3 (r )S p (r )dr Q - -^^ r ^ +1 s a (l.O)S (l.O)) 
Jo £> r o 27+I 

+A f G s (r,r )S 7 {r )dr , (95) 
Jo 

T 7 (r) = y^Rmif G t (r,r )F 4 {r )Tp(r )dr + f G t (r,r )F 5 (r )Sp(r )dr 

a/3 J ° J ° 

1 dGtir > r0) F 6 (ro)S p (ro)dr + t ^pAp^T^dr, 



with 



dr a Jo dr 

-c 10 rT+ 1 SQ (1.0)S'/,(1.0) + ^5 Q (r)5 / 3(r)]+A / Gt^^T^dro, (96) 

r ./n 



7-n d { Sat \ . 1 ds a 

^1 = Ci— - c 3 — ( — ) + c 4 — — , 
To dr dr Q 

TP Sa 

r-X = C2-2, 

^3 = -C3-2, 
^0 

t Q d , s a . 1 . c?s Q 2 . 1 ds a 

= - C 8 — (-J ) + C8jh Sa)+Cg-2- , 

dr rg rg dr r rg dr 

F5 = - C6 ^ ( | ) + C7 4 ( ^"^ q) + C107(7 + 1) ^ + C10 ^! ( S ) 



d ,1 ds (y » d , §qi \ "\. , d s 2 ds^ 

13— (~2 "J— J - C12 3— (-3) + c 13-2-(-7^ 3— 

£ a (i , S a . 1 ds a 5 a 
-9 + ■'Cio— — 57 — Cii-g— C\2—^i 

rg dr rg rg dr rg 



j-, Sa 
^7 = - c 8 — 
r 



The discretization of this integral equation system is done along the lines described in subsection IV Al 
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TABLE IV: Convergence of the integral equation approach (IEA) and the Dudley and James method (D&J). We show the 
dependence of A on the radial grid number N for Rm = 50 and Rm = 80 using a truncation degree L — 9. The interpolation 
in the IEA case is done with a fit of the date to a function a + bN~ 2 . Dudley and James had used Richardson extrapolation 
based on the values for N = 75, 100, 125. 







N=50 


N=75 


N=100 


N=125 


Extrapolation 


IEA Rm = 


50 


-23.50+4.97i 


-22.97+4.03i 


-22.76+3.66i 


-22.63+3.34i 




IEA Rm = 


SO 


-29.35+9.46i 


-27.70+7.72i 


26.95+6.98i 


-26.51+6.57i 


-26.10+ 6.26i 


D&J Rm = 


80 




26.32+6.01i 


26.33+6.04i 


26.33+6.05i 


-26.34+6.06i 



B. A numerical example - The Bullard-Gellman model 

In the following, we will test the suitability of the integral equation approach to the simulation of large scale velocity 
fields for a particular dynamo model. In 1954, Bullard and Gellman |31| had studied the flow structure 

u = sf + 5t° (97) 

where 

sl c {r)=r\l-rf (98) 

f°(r) =r 2 (l-r) (99) 

claiming that this flow acts indeed as a dynamo. Later, using higher spatial resolution, Gibson and Roberts |32j | and 
Dudley and James [33] falsified this result showing that there is no dynamo up to a magnetic Reynolds number of 80. 

Here we treat the Bullard-Gellman model within the framework of the integral equation approach. 

In Fig. 0we plot the real and the imaginary part of the first eigenvalue of the Dl solution (in the terminology of 
Dudley and James) for the Bullard-Gellman dynamo model in dependence on the magnetic Reynolds number. The 
truncation degree is L — 9, the number of radial grid points is N — 75. This is essentially the same curve as published 
in [3^ where a truncation degree L = 12 and a number of grid points of N = 100 had been used, however. 

In Table (|IV|I we have compiled some results concerning the convergence of our method and the method of Dudley 
and James. For Rm = 50 (where no data are available from Dudley and James) we see a reasonable convergence of 
the real part but a slow convergence of the imaginary part. The latter might be due to the fact that we are not far 
from the transition point to oscillatory behaviour where the imaginary part is sensible to changes in the grid number. 
For Rm = 80 we have to concede that the convergence in our case is slower than in the differential equation method 
of Dudley and James. 

A similar conclusion can be drawn from the treatment of other models (Lilley model, modified Lilley model). 
Although our method yields essentially the same results as the differential equation approach, it seems worth to look 
for refined numerical methods to solve the integral eigenvalue equation. 

VII. REMARKS AND CONCLUSIONS 

We have established the integral equation approach to time-dependent kinematic dynamos in arbitrary domains. 
This approach is based on Biot-Savart's law. The main advantage of the method is its suitability to handle dynamos in 
arbitrary domains. The necessity to solve the Laplace equation in the exteriour of the dynamo domain is circumvented 
by the (implicit) solution of boundary integral equations for the electric potential and the magnetic vector potential. 

It should be noted that we have worked out only one possible form of the integral equation approach which results in 
a linear eigenvalue problem. Another form could start from rewriting Eq. I|17fl into the form of a Helmholtz equation 
for the vector potential 

AA - n a\A = /i er(F - V<fi) . (100) 
Then, the pendant to Eq. i|19[l would read 

r'|) dV 

-fc|r-r'|) dS', (101) 



F(r') x (r - r') 



D 



r — r 



/|3 



-— / V (B)n(B)x 



exp (k\r — r'|)(l — fc|r - 
exp(fc|r-r'|)(l 



r — s 
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FIG. 6: Real and imaginary part of the first eigenvalue for the Bullard-Gellman dynamo model in dependence on the magnetic 
Reynolds number. The truncation degree is L — 9, the number of radial grid points is N = 75. 



with 



(102) 



Without going into the details of such a formulation (we skip here the equations for the electric potential and the 
vector potential at the boundary), we see immediately that we end up with a non- linear eigenvalue equation for the 
eigenvalue A. It would be interesting to compare the numerical performance of such a formulation with the present 
one. 

We plan to use our formulation for a number of dynamo problems, in particular problems which are connected with 
the design and optimization for new dynamo experiments and with velocity reconstruction problem for the existing 
experiments. 



VIII. APPENDIX A 



In this part, we give another approach to establish the integral equations l|47|) and (|67|l . We start with the following 
differential equation problem: 



Xsi 



1 Am _ 1(1 + 1) 

/ioc dr 2 r 2 



Sim] + a{r)tir, 



1 ,d 2 ti m 1(1 + 1) d . ,dsi m 1(1 + 1) 

Xtlm = — 7^ o Urn - ~r( a ( r )— j — H 5 Ot[r)si m , 

fio<j ar z r z dr dr r l 

t lm (R) = R^r L \ r = R + ls lm (R)=0, 
dr 



(103) 
(104) 
(105) 



First, we derive the Green's function, G s (r, ro), corresponding to the equation H103|l . This Green's function satisfies: 

(106) 
(107) 
(108) 



d 2 G s 
dr 2 


Hl + l) r 

r 2 s 


= S(r - r 




G s \ r —Q 


= o, 


R dGs \ 
R dr lr = 


=R + lG s \ r= R 


= . 
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According to the construction method of Green's functions ([3(|, P- 355), we obtain G s (r, fo) in the following form: 

g ^o) = \ zw^y^ dog) 

[ 2/+1 'o ' ' ' - 'o- 

As for Eq. (|104fl . we first rewrite it in the form: 

1 d 2 F 1 1(1 + 1) d Aa s , . 

[Iqo dr z [i$o r A dr dr 

where F = ti m — fiQaasi m . This differential equation problem for ti m can be split into two problems. One of them 
reads 

d 2 F 1 l(r+l) d da 

-7-5 5 Fl + f-lO<T—(S[ m —) - ^ a\ti rn = 0, (111) 

dr z r z dr dr 

Fl\ r =R = 0. 



The other is 

d 2 F 2 1(1 + 1) 



F 2 = 0, (112) 



dr 2 r 
F 2 \ r =B. = -Hocra(R)si m (R). 

For the differential equation problem (|lllfl . applying the construction method of Green's function (|30f) again, we 
obtain its Green's function in the form: 

Gt(r,r ) = { ^j^-Sp;^' (113) 

which satisfies 

G*|r=o = 0, (115) 

G t \r=R = 0. (116) 

As for the differential equation problem l|112|) . the solution can be expressed as: 

F2 = -j^a(R)s lm (R)r l+ \ (117) 

Then the superposition theorem of the linear problems allows us to obtain the following integral equations for s; m 
and t hn : 

pR rR. 

sim(r) = G s (r,ro)noaa(ro)ti m (ro)dro + noaX Gs(r,r )si m (r )dro, (118) 

Jo Jo 

tim(r) = fi a-a(r)s lm (r) - / ii a^(si m (r )^r^-)Gt(r, r )dr + 

Jo dr o dr 

i-R r i+i 
Xfi a J G t (r,r a )t lm (r )dr - -j^ii aa(R)s lm (R). (119) 

Integrating by parts the terms containing derivatives of s; m in Eq. 1|119|) . we obtain 

rR rR. 

si m (r) = - G s (r, r )iM)cra(ro)ti m (r )dro + ^ofA / G s {r,r )s hn (r )dro, (120) 



, i s , , ( \ \ f da{r ) dG t {r,r ) 

tim{r) = (ioaa{r)s[ m (r) + / //o^— s lm (r ) dr 

Jo dr o or 

f R r l+1 
+A M oa J Gt(r,r a )ti m (r a )dr a - ^^■n aa(R)si m (R). (121) 

Therefore, we have obtained the same integral equations as expressed in Eqs. (|47|l and Q6 7[). 
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